Fit models with metabolic index and temperature velocities for each species
Max Lindmark
20 March, 2023
Intro
Main analysis: This scripts reads spatially referenced linear trends over time in biomass, metabolic index and temperature (“03_calculate_velocities.Rmd”), and fits spatial sdmTMB models to those
# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(tidylog)
library(RCurl)
library(devtools)
library(sdmTMB)
library(rMR)
library(patchwork)
library(raster)
library(VoCC)
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(viridis)
library(ggh4x)
library(ggridges)
# Source code for plots
source_url("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/R/functions/map-plot.R")
#> Warning: package 'sf' was built under R version 4.0.5
theme_set(theme_plot())
pal <- brewer.pal(name = "Dark2", n = 8)Read and plot data
# Read from GH?
d <- read_csv("data/clean/velocities.csv") %>%
mutate(quarter_f = as.factor(quarter))
#> Rows: 87885 Columns: 20
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (4): id, species, life_stage, group
#> dbl (16): X, Y, quarter, mean_log_biomass, mean_temp, mean_oxy, mean_phi, me...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: new variable 'quarter_f' (factor) with 2 unique values and 0% NA
glimpse(d)
#> Rows: 87,885
#> Columns: 21
#> $ id <chr> "320.016716551297_6028.58668110674_cod_adult_1", "…
#> $ X <dbl> 320.0167, 320.0167, 320.0167, 320.0167, 320.0167, …
#> $ Y <dbl> 6028.587, 6028.587, 6028.587, 6028.587, 6028.587, …
#> $ species <chr> "cod", "cod", "cod", "cod", "dab", "dab", "dab", "…
#> $ life_stage <chr> "adult", "adult", "juvenile", "juvenile", "adult",…
#> $ quarter <dbl> 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1,…
#> $ mean_log_biomass <dbl> 4.3916628, 4.5915365, 2.8446881, 2.9603249, 3.6929…
#> $ mean_temp <dbl> 3.369871, 10.260329, 3.369871, 10.260329, 3.369871…
#> $ mean_oxy <dbl> 8.114926, 6.209476, 8.114926, 6.209476, 8.114926, …
#> $ mean_phi <dbl> 1.5023498, 0.8518787, 4.6064940, 2.6120242, 4.4595…
#> $ mean_log_biomass_ct <dbl> 1.545992832, 1.745866571, -0.000981832, 0.11465494…
#> $ mean_temp_ct <dbl> -2.888372, 4.002085, -2.888372, 4.002085, -2.88837…
#> $ mean_oxy_ct <dbl> 2.3313114, 0.4258614, 2.3313114, 0.4258614, 2.3313…
#> $ mean_phi_ct <dbl> -1.6917880, -2.3422591, 1.4123562, -0.5821135, 1.2…
#> $ bio_vel <dbl> -0.03357378, -0.03357378, -0.03601321, -0.03601321…
#> $ phi_vel <dbl> 0.3286411, 0.5693023, 0.3286411, 0.5693023, 0.3286…
#> $ temp_vel <dbl> 1.6318561, 0.7592320, 1.6318561, 0.7592320, 1.6318…
#> $ group <chr> "cod_adult", "cod_adult", "cod_juvenile", "cod_juv…
#> $ phi_vel_sc <dbl> 0.8134520, 0.9991260, 0.8942800, 1.1221752, 1.0096…
#> $ temp_vel_sc <dbl> -0.23950578, -0.47555486, -0.11759968, -0.39684290…
#> $ quarter_f <fct> 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1,…
# Biotic velocities
plot_map +
geom_raster(data = d, aes(X*1000, Y*1000, fill = bio_vel)) +
facet_grid(life_stage ~ species) +
#scale_fill_viridis() +
scale_fill_gradient2() +
geom_sf(size = 0.3)
#> Warning: Removed 2425 rows containing missing values (geom_raster).
# Phi velocities
plot_map +
geom_raster(data = d, aes(X*1000, Y*1000, fill = phi_vel)) +
facet_grid(life_stage ~ species) +
#scale_fill_viridis() +
scale_fill_gradient2() +
geom_sf(size = 0.3)
#> Warning: Removed 2425 rows containing missing values (geom_raster).
# Temperature velocities
plot_map +
geom_raster(data = d, aes(X*1000, Y*1000, fill = temp_vel)) +
facet_grid(life_stage ~ species) +
#scale_fill_viridis() +
scale_fill_gradient2() +
geom_sf(size = 0.3)
#> Warning: Removed 2425 rows containing missing values (geom_raster).
# Plot biotic vs phi and temperature velocities
ggplot(d, aes(phi_vel, bio_vel)) +
geom_point(alpha = 0.4) +
facet_grid(life_stage ~ species) +
stat_smooth(method = "lm") +
theme(aspect.ratio = 1) +
NULL
#> `geom_smooth()` using formula 'y ~ x'
ggplot(d, aes(temp_vel, bio_vel)) +
geom_point(alpha = 0.4) +
facet_grid(life_stage ~ species) +
stat_smooth(method = "lm") +
theme(aspect.ratio = 1) +
NULL
#> `geom_smooth()` using formula 'y ~ x'
# Only cod
head(d) %>% as.data.frame()
#> id X Y species
#> 1 320.016716551297_6028.58668110674_cod_adult_1 320.0167 6028.587 cod
#> 2 320.016716551297_6028.58668110674_cod_adult_4 320.0167 6028.587 cod
#> 3 320.016716551297_6028.58668110674_cod_juvenile_1 320.0167 6028.587 cod
#> 4 320.016716551297_6028.58668110674_cod_juvenile_4 320.0167 6028.587 cod
#> 5 320.016716551297_6028.58668110674_dab_adult_1 320.0167 6028.587 dab
#> 6 320.016716551297_6028.58668110674_dab_adult_4 320.0167 6028.587 dab
#> life_stage quarter mean_log_biomass mean_temp mean_oxy mean_phi
#> 1 adult 1 4.391663 3.369871 8.114926 1.5023498
#> 2 adult 4 4.591537 10.260329 6.209476 0.8518787
#> 3 juvenile 1 2.844688 3.369871 8.114926 4.6064940
#> 4 juvenile 4 2.960325 10.260329 6.209476 2.6120242
#> 5 adult 1 3.692946 3.369871 8.114926 4.4595684
#> 6 adult 4 3.923447 10.260329 6.209476 2.5287129
#> mean_log_biomass_ct mean_temp_ct mean_oxy_ct mean_phi_ct bio_vel
#> 1 1.545992832 -2.888372 2.3313114 -1.6917880 -0.03357378
#> 2 1.745866571 4.002085 0.4258614 -2.3422591 -0.03357378
#> 3 -0.000981832 -2.888372 2.3313114 1.4123562 -0.03601321
#> 4 0.114654948 4.002085 0.4258614 -0.5821135 -0.03601321
#> 5 0.847275563 -2.888372 2.3313114 1.2654306 0.21900577
#> 6 1.077776951 4.002085 0.4258614 -0.6654249 0.21900577
#> phi_vel temp_vel group phi_vel_sc temp_vel_sc quarter_f
#> 1 0.3286411 1.631856 cod_adult 0.813452 -0.2395058 1
#> 2 0.5693023 0.759232 cod_adult 0.999126 -0.4755549 4
#> 3 0.3286411 1.631856 cod_juvenile 0.894280 -0.1175997 1
#> 4 0.5693023 0.759232 cod_juvenile 1.122175 -0.3968429 4
#> 5 0.3286411 1.631856 dab_adult 1.009634 -0.2510603 1
#> 6 0.5693023 0.759232 dab_adult 1.335121 -0.5594185 4
d_cod <- d %>%
filter(group == "cod_adult")
#> filter: removed 70,544 rows (80%), 17,341 rows remaining
p1 <- plot_map +
geom_raster(data = d_cod, aes(X*1000, Y*1000, fill = bio_vel)) +
scale_fill_gradient2(name = "Biotic velocity") +
geom_sf(size = 0.3)
p2 <- plot_map +
geom_raster(data = d_cod, aes(X*1000, Y*1000, fill = temp_vel)) +
scale_fill_gradient2(name = "Temperature velocity") +
geom_sf(size = 0.3)
p3 <- plot_map +
geom_raster(data = d_cod, aes(X*1000, Y*1000, fill = phi_vel)) +
scale_fill_gradient2(name = "\u03C6 velocity") +
geom_sf(size = 0.3)
((p2 + p3) / (p1 + plot_spacer()))
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.
ggsave("figures/map_velos_cod.pdf", width = 17, height = 17, units = "cm", device = cairo_pdf)
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted. Consider using geom_tile() instead.
#> Raster pixels are placed at uneven vertical intervals and will be shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.Fit sdmTMB models of biomass velocities as a function of phi velocities
coefs <- list()
pred_velo1 <- list()
pred_velo2 <- list()
for(i in unique(d$group)){
dd <- filter(d, group == i)
mesh <- make_mesh(dd, c("X", "Y"), n_knots = 200)
plot(mesh)
m <- sdmTMB(
bio_vel ~ phi_vel_sc * mean_phi_ct + mean_log_biomass_ct + quarter_f,
mesh = mesh,
data = dd,
family = gaussian(),
spatial = "on",
spatiotemporal = "off",
control = sdmTMBcontrol(newton_loops = 1)
)
sanity(m)
# Save coefficients??
coefs[[i]] <- tidy(m, conf.int = TRUE) %>% mutate(group = i)
# Now predict the biotic velocity for high and low temperature with a sd change in the climate velocity
low_mean_clim <- as.numeric(quantile(dd$mean_phi_ct, prob = 0.05))
high_mean_clim <- as.numeric(quantile(dd$mean_phi_ct, prob = 0.95))
#sd_clim_v <- sd(dd$phi_vel)
mean_log_biomass_ct <- mean(dd$mean_log_biomass_ct)
# Predict
nd <- data.frame(phi_vel_sc = c(-1, 0, -1, 0),
mean_phi_ct = rep(c(low_mean_clim, high_mean_clim), each = 2),
mean_log_biomass_ct = rep(mean_log_biomass_ct, 4),
X = mean(dd$X),
Y = mean(dd$Y),
quarter_f = as.factor(4))
nsim = 1000
pred <- predict(m, newdata = nd, se_fit = TRUE, nsim = nsim)
pred <- data.frame(
mean_phi_ct = rep(c(low_mean_clim, high_mean_clim), each = nsim),
est_change = c(pred[1, ], pred[3, ]),
est_0 = c(pred[2, ], pred[4, ]))
pred <- pred %>%
mutate(bio_trend_delta = est_change - est_0, # if bio trend is positive, bio increases with declines in phi
mean_clim = ifelse(mean_phi_ct == low_mean_clim,
"low_phi", "high_phi"),
group = i)
pred_velo1[[i]] <- pred
# Now make a prediction across a range of velocities
# Again I'm predicting for the range of the specific model
nd3 <- data.frame(expand.grid(phi_vel_sc = seq(quantile(dd$phi_vel_sc, probs = 0.05),
quantile(dd$phi_vel_sc, probs = 0.95),
length.out = 50),
mean_phi_ct = c(low_mean_clim, high_mean_clim))) %>%
mutate(mean_log_biomass_ct = mean_log_biomass_ct,
X = mean(dd$X),
Y = mean(dd$Y),
quarter_f = as.factor(4))
pred3 <- predict(m, newdata = nd3, se_fit = TRUE)
pred3 <- pred3 %>%
mutate(lwr = est - 1.96*est_se,
upr = est + 1.96*est_se,
group = i)
pred_velo2[[i]] <- pred3
ggplot(pred3, aes(phi_vel_sc, est, color = factor(mean_phi_ct))) + geom_line()
}
#> filter: removed 70,544 rows (80%), 17,341 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 72,173 rows (82%), 15,712 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 86,033 rows (98%), 1,852 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 86,343 rows (98%), 1,542 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 64,646 rows (74%), 23,239 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 79,901 rows (91%), 7,984 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 79,755 rows (91%), 8,130 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 75,800 rows (86%), 12,085 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#> new variable 'upr' (double) with 100 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
phi_coefs <- dplyr::bind_rows(coefs) %>%
mutate(sig = "N",
sig = ifelse(estimate < 0 & conf.high < 0, "Y", sig),
sig = ifelse(estimate > 0 & conf.low > 0, "Y", sig))
#> mutate: new variable 'sig' (character) with 2 unique values and 0% NA
phi_pred_delta <- dplyr::bind_rows(pred_velo1) %>%
separate(group, c("species", "life_stage"),
sep = "_", remove = FALSE) %>%
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage))
#> mutate: changed 16,000 values (100%) of 'species' (0 new NA)
#> changed 16,000 values (100%) of 'life_stage' (0 new NA)
phi_pred <- dplyr::bind_rows(pred_velo2) %>%
separate(group, c("species", "life_stage"),
sep = "_", remove = FALSE) %>%
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage))
#> mutate: changed 800 values (100%) of 'species' (0 new NA)
#> changed 800 values (100%) of 'life_stage' (0 new NA)
plot_map +
geom_raster(data = d, aes(X*1000, Y*1000, fill = bio_vel)) +
facet_grid(life_stage ~ species) +
#scale_fill_viridis() +
scale_fill_gradient2() +
geom_sf(size = 0.3)
#> Warning: Removed 2425 rows containing missing values (geom_raster).
write_csv(phi_coefs, "output/phi_velo_coefs.csv")
write_csv(phi_pred_delta, "output/pred_phi_velo_delta.csv")
write_csv(phi_pred, "output/pred_phi_velo.csv")Plot phi velocity models
phi_pred_delta2 <- phi_pred_delta %>%
mutate(mean_clim = ifelse(mean_clim == "high_phi", "High \u03C6", "Low \u03C6")) %>%
mutate(id2 = paste(species, life_stage))
#> mutate: changed 16,000 values (100%) of 'mean_clim' (0 new NA)
#> mutate: new variable 'id2' (character) with 8 unique values and 0% NA
order <- phi_pred_delta2 %>%
group_by(id2, mean_clim) %>%
summarise(mean_delta = mean(bio_trend_delta)) %>%
pivot_wider(names_from = mean_clim, values_from = mean_delta) %>%
mutate(diff = `Low φ` - `High φ`) %>%
arrange(desc(diff)) %>%
mutate(Expected = ifelse(`Low φ` < `High φ`, "Yes", "No"))
#> group_by: 2 grouping variables (id2, mean_clim)
#> summarise: now 16 rows and 3 columns, one group variable remaining (id2)
#> pivot_wider: reorganized (mean_clim, mean_delta) into (High φ, Low φ) [was 16x3, now 8x3]
#> mutate (grouped): new variable 'diff' (double) with 8 unique values and 0% NA
#> mutate (grouped): new variable 'Expected' (character) with 2 unique values and 0% NA
order
#> # A tibble: 8 × 5
#> # Groups: id2 [8]
#> id2 `High φ` `Low φ` diff Expected
#> <chr> <dbl> <dbl> <dbl> <chr>
#> 1 Dab Juvenile 0.0154 0.0929 0.0775 No
#> 2 Dab Adult 0.0252 0.0382 0.0130 No
#> 3 Cod Juvenile 0.0226 0.0335 0.0108 No
#> 4 Cod Adult 0.00753 -0.0102 -0.0177 Yes
#> 5 Flounder Adult 0.0293 -0.0118 -0.0412 Yes
#> 6 Flounder Juvenile 0.101 -0.123 -0.224 Yes
#> 7 Plaice Juvenile 0.141 -0.0893 -0.230 Yes
#> 8 Plaice Adult 0.0700 -0.178 -0.247 Yes
phi_pred_delta2 <- phi_pred_delta2 %>%
left_join(order %>% dplyr::select(Expected), by = "id2")
#> Adding missing grouping variables: `id2`
#> left_join: added one column (Expected)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 16,000
#> > ========
#> > rows total 16,000
phi_pred_delta2$id2 <- factor(phi_pred_delta2$id2, levels = c(order$id2))
phi_pred_delta2 %>%
ggplot(aes(bio_trend_delta, id2, fill = mean_clim, alpha = Expected)) +
geom_vline(xintercept = 0, linetype = 2, color = "grey30") +
geom_density_ridges(color = NA) +
scale_fill_manual(values = c("steelblue2", "tomato3")) +
scale_alpha_manual(values = c(0.4, 0.8)) +
labs(y = "", x = "\u0394 biotic velocity with a 1 st.dev\ndecline in \u03C6 velocity", fill = "Mean climate") +
theme_plot(base_size = 16) +
theme(legend.direction = "vertical") +
NULL
#> Picking joint bandwidth of 0.0071
ggsave("figures/mi_velo_delta.pdf", width = 17, height = 17, units = "cm", device = cairo_pdf)
#> Picking joint bandwidth of 0.0071
# phi_pred %>%
# mutate(id2 = paste(life_stage, species)) %>%
# group_by(group) %>%
# mutate(mean_clim = ifelse(mean_phi_ct == max(mean_phi_ct), "High \u03C6", "Low \u03C6")) %>%
# ungroup() %>%
# ggplot(aes(phi_vel_sc, est, color = mean_clim, fill = mean_clim)) +
# geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.3, color = NA) +
# geom_line(alpha = 0.8) +
# labs(y = "Predicted Biotic velocity (km/year)",
# x = "\u03C6 velocity (km/year)") +
# scale_fill_manual(values = c("steelblue2", "tomato3"), name = "") +
# scale_color_manual(values = c("steelblue2", "tomato3"), name = "") +
# facet_grid2(life_stage ~ species, scales = "free_y", independent = "y") +
# theme(aspect.ratio = 1) +
# NULL
#ggsave("figures/mi_velo_conditional.pdf", width = 17, height = 17, units = "cm", device = cairo_pdf)Fit sdmTMB models of biomass velocities as a function of temperature velocities, given mean temperature (ooxygen would be cool to but I think perhaps its too correlated with temperature, cor = 0.8 ish)
pred_velo_temp <- list()
pred_velo_temp2 <- list()
for(i in unique(d$group)){
dd <- filter(d, group == i)
mesh <- make_mesh(dd, c("X", "Y"), n_knots = 200)
plot(mesh)
m <- sdmTMB(
bio_vel ~ temp_vel_sc * mean_temp_ct + mean_log_biomass_ct + quarter_f,
mesh = mesh,
data = dd,
family = gaussian(),
spatial = "on",
spatiotemporal = "off",
control = sdmTMBcontrol(newton_loops = 1)
)
sanity(m)
# Predict
low_mean_clim <- as.numeric(quantile(dd$mean_temp_ct, prob = 0.05))
high_mean_clim <- as.numeric(quantile(dd$mean_temp_ct, prob = 0.95))
mean_log_biomass_ct <- mean(dd$mean_log_biomass_ct)
nd <- data.frame(expand.grid(temp_vel_sc = c(0, 1),
mean_temp_ct = seq(low_mean_clim,
high_mean_clim,
length.out = 2))) %>%
mutate(mean_log_biomass_ct = mean_log_biomass_ct,
X = mean(dd$X),
Y = mean(dd$Y),
quarter_f = as.factor(4))
nsim <- 1000
pred <- predict(m, newdata = nd, nsim = nsim)
pred <- data.frame(
mean_temp_ct = rep(c(low_mean_clim, high_mean_clim), each = nsim),
est_0 = c(pred[1, ], pred[3, ]),
est_change = c(pred[2, ], pred[4, ])) # note different order here than phi
pred <- pred %>%
mutate(delta_bio_vel = est_change - est_0, # if bio trend is positive, bio increases with declines in phi
mean_clim = ifelse(mean_temp_ct == low_mean_clim,
"low_temperature", "high_temperature"),
group = i)
pred_velo_temp[[i]] <- pred
# Predict 2 (max and min warming, what's the expected?)
# nd2 <- data.frame(expand.grid(temp_vel_sc = c(quantile(dd$temp_vel_sc, probs = 0.01),
# quantile(dd$temp_vel_sc, probs = 0.99)),
# mean_temp_ct = seq(quantile(dd$mean_temp_ct, probs = 0.01),
# quantile(dd$mean_temp_ct, probs = 0.99),
# length.out = 2))) %>%
# mutate(mean_log_biomass_ct = 0,
# X = mean(dd$X),
# Y = mean(dd$Y),
# quarter_f = as.factor(4))
#
# pred2 <- predict(m, newdata = nd2) %>% mutate(group = i)
#
# pred_velo_temp2[[i]] <- pred2
}
#> filter: removed 70,544 rows (80%), 17,341 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'delta_bio_vel' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 72,173 rows (82%), 15,712 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'delta_bio_vel' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 86,033 rows (98%), 1,852 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'delta_bio_vel' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 86,343 rows (98%), 1,542 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'delta_bio_vel' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 64,646 rows (74%), 23,239 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'delta_bio_vel' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 79,901 rows (91%), 7,984 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✖ `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> ✖ `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> ✖ `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'delta_bio_vel' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 79,755 rows (91%), 8,130 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'delta_bio_vel' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 75,800 rows (86%), 12,085 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#> new variable 'X' (double) with one unique value and 0% NA
#> new variable 'Y' (double) with one unique value and 0% NA
#> new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'delta_bio_vel' (double) with 2,000 unique values and 0% NA
#> new variable 'mean_clim' (character) with 2 unique values and 0% NA
#> new variable 'group' (character) with one unique value and 0% NA
pred_velo_temp <- dplyr::bind_rows(pred_velo_temp) %>%
separate(group, c("species", "life_stage"),
sep = "_", remove = FALSE) %>%
mutate(species = str_to_title(species),
life_stage = str_to_title(life_stage))
#> mutate: changed 16,000 values (100%) of 'species' (0 new NA)
#> changed 16,000 values (100%) of 'life_stage' (0 new NA)
# pred_velo_temp2 <- dplyr::bind_rows(pred_velo_temp2) %>%
# separate(group, c("species", "life_stage"),
# sep = "_", remove = FALSE) %>%
# mutate(species = str_to_title(species),
# life_stage = str_to_title(life_stage))
write_csv(pred_velo_temp, "output/pred_temp_velo.csv")pred_velo_temp2 <- pred_velo_temp %>%
mutate(mean_clim = ifelse(mean_clim == "low_temperature", "Cold", "Warm")) %>%
mutate(id2 = paste(species, life_stage))
#> mutate: changed 16,000 values (100%) of 'mean_clim' (0 new NA)
#> mutate: new variable 'id2' (character) with 8 unique values and 0% NA
order <- pred_velo_temp2 %>%
group_by(id2, mean_clim) %>%
summarise(mean_delta = mean(delta_bio_vel)) %>%
pivot_wider(names_from = mean_clim, values_from = mean_delta) %>%
mutate(diff = `Warm` - `Cold`) %>%
arrange(desc(diff)) %>%
mutate(Expected = ifelse(Cold > Warm, "Yes", "No"))
#> group_by: 2 grouping variables (id2, mean_clim)
#> summarise: now 16 rows and 3 columns, one group variable remaining (id2)
#> pivot_wider: reorganized (mean_clim, mean_delta) into (Cold, Warm) [was 16x3, now 8x3]
#> mutate (grouped): new variable 'diff' (double) with 8 unique values and 0% NA
#> mutate (grouped): new variable 'Expected' (character) with 2 unique values and 0% NA
order
#> # A tibble: 8 × 5
#> # Groups: id2 [8]
#> id2 Cold Warm diff Expected
#> <chr> <dbl> <dbl> <dbl> <chr>
#> 1 Flounder Adult -0.0226 0.106 0.128 No
#> 2 Dab Adult 0.0289 0.0553 0.0264 No
#> 3 Dab Juvenile 0.0142 0.0313 0.0171 No
#> 4 Plaice Juvenile 0.0424 0.0509 0.00852 No
#> 5 Plaice Adult -0.0250 -0.0382 -0.0132 Yes
#> 6 Cod Adult 0.0330 -0.00467 -0.0377 Yes
#> 7 Flounder Juvenile 0.0488 -0.0242 -0.0730 Yes
#> 8 Cod Juvenile 0.0410 -0.0358 -0.0769 Yes
pred_velo_temp2 <- pred_velo_temp2 %>%
left_join(order %>% dplyr::select(Expected), by = "id2")
#> Adding missing grouping variables: `id2`
#> left_join: added one column (Expected)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 16,000
#> > ========
#> > rows total 16,000
pred_velo_temp2$id2 <- factor(pred_velo_temp2$id2, levels = c(order$id2))
pred_velo_temp2 %>%
ggplot(aes(delta_bio_vel, id2, fill = mean_clim, alpha = Expected)) +
geom_vline(xintercept = 0, linetype = 2, color = "grey30") +
geom_density_ridges(color = NA) +
scale_fill_manual(values = c("steelblue2", "tomato3")) +
scale_alpha_manual(values = c(0.4, 0.8)) +
labs(y = "", x = "\u0394 biotic velocity with a 1 st.dev\nincrease in temperature velocity", fill = "Mean temperature") +
theme_plot(base_size = 16) +
theme(legend.direction = "vertical") +
NULL
#> Picking joint bandwidth of 0.00598
ggsave("figures/temp_velo_delta.pdf", width = 15, height = 17, units = "cm", device = cairo_pdf)
#> Picking joint bandwidth of 0.00598knitr::knit_exit
#> function (append, fully = TRUE)
#> {
#> if (missing(append))
#> append = if (out_format(c("latex", "sweave", "listings")))
#> "\\end{document}"
#> else if (out_format("html"))
#> "</body>\n</html>"
#> else ""
#> .knitEnv$terminate = append
#> .knitEnv$terminate_fully = fully
#> invisible()
#> }
#> <bytecode: 0x7fa2ad793708>
#> <environment: namespace:knitr># pred_velo_temp2 %>%
# mutate(group2 = paste(species, life_stage, sep = "")) %>%
# group_by(group) %>%
# mutate(mean_clim = ifelse(mean_temp_ct == min(mean_temp_ct), "Low temperature", "High temperature"),
# warming = ifelse(temp_vel_sc == max(temp_vel_sc), "Max warming", "Min warming")) %>%
# ungroup() %>%
# ggplot(aes(est, group, color = mean_clim)) +
# geom_point(alpha = 0.7, size = 3) +
# scale_color_manual(values = c("red", "blue")) +
# geom_vline(xintercept = 0, linetype = 2) +
# facet_wrap(~factor(warming, levels = c("Min warming", "Max warming"))) +
# scale_fill_gradient2()
#
# # Filter only max warming
# pred_velo_temp2 %>%
# mutate(group2 = paste(species, life_stage, sep = "")) %>%
# group_by(group) %>%
# mutate(mean_clim = ifelse(mean_temp_ct == min(mean_temp_ct), "Low temperature", "High temperature"),
# warming = ifelse(temp_vel_sc == max(temp_vel_sc), "Max warming", "Min warming")) %>%
# ungroup() %>%
# filter(warming == "Max warming") %>%
# ggplot(aes(est, group, color = mean_clim)) +
# geom_point(alpha = 0.7, size = 3) +
# scale_color_manual(values = c("red", "blue")) +
# geom_vline(xintercept = 0, linetype = 2) +
# facet_wrap(~factor(warming, levels = c("Min warming", "Max warming"))) +
# scale_fill_gradient2()